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Figure  1.  Bony  anatomy  of  the  spinal  column  (panel  a)  and  a  typical  vertebra  (panel  b). 

Vertebra  are  color-coded  according  to  their  location  classification.  Panel  c  is  an 
illustration  (not  drawn  to  scale)  of  an  intervertebral  disc  showing  the  lamellar 
architecture  of  the  annulus  fibrosus  (white)  which  surrounds  the  nucleus  pulposus 
(grey).  Some  layers  have  been  cut  away  from  the  annulus  fibrosus  to  show  the  fiber 
network  within  the  lamellae.  Note  that  only  four  layers  of  lamellae  are  depicted  in  the 
figure  but  the  annulus  fibrosus  usually  has  15-25  layers.  Collagen  fibers  (diagonal 
lines)  are  oriented  at  ±30°  to  the  transverse  plane  of  the  disc,  with  the  direction 
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Tzz  (black),  Txy  (cyan),  Ty^  (magenta),  and  (green)  for  three  sets  of  fiber 
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orientation  vectors  ao.  Panel  a  shows  the  stress  response  of  a  single  fiber  family  with 
initial  direction  vector  in  the  Y Z-plane.  Similarly,  panels  b  and  c  show  a  fiber  family 
with  initial  direction  vector  in  the  XZ-plane,  and  XY -plane,  respectively . 14 


V 


Figure  5.  Two  fiber  families  compression  test.  Simulation  (symbols)  comparison  against 
theoretical  (solid  lines)  for  the  components  of  the  Cauchy  stress  T^x  (rcd),  Tyy  (blue), 

Tzz  (black),  Txy  (cyan),  Tyz  (magenta),  and  T^x  (green)  for  two  liber  family 

orientation  vectors  and  Qq .  16 

Figure  6.  Local  coordinate  system  for  an  intervertebral  disc.  An  idealized  intervertebral 
disc  in  the  reference  coordinate  system.  The  point  r  is  on  the  surface  of  the  cylinder 
with  corresponding  normal  vector  n,  chosen  tangent  vector  t,  and  binormal  vector  b. 

The  vectors  ao  and  Qq  for  this  point  are  also  shown .  17 
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1.  Introduction 


There  is  an  increasing  interest  within  the  automotive,  military,  and  clinical  communities  to  use 
computer  simulations  to  predict  injury  in  humans.  Researchers  have  developed  finite  element 
models  of  biological  tissues  and  subjected  their  models  to  various  loading  conditions  (1^). 

Finite  element  models  rely  on  accurate  constitutive  models  of  these  tissues  and  therefore 
researchers  must  make  critical  decisions  regarding  the  appropriate  level  of  detail  to  include. 
While  it  is  safe  to  assume  anisotropy  does  exist  in  biological  tissues,  the  role  it  plays  during  a 
high  loading  rate  event,  such  as  a  blast,  is  not  well  understood.  Numerical  models  that  can  turn 
these  effects  on  or  off  are  an  incredibly  useful  tool. 

The  biological  tissues  that  we  consider  here  are  soft  and  often  have  a  high  water  content  that 
places  their  bulk  moduli  close  to  that  of  water,  i.e.,  around  2.3  GPa.  However,  they  also  typically 
have  relatively  low  shear  moduli,  making  them  nearly  incompressible.  An  extreme  example  is 
brain  tissue  where  the  lowest  shear  modulus  values  measured  are  around  2  kPa,  cf.  (5-10). 
Because  these  tissues  have  low  shear  moduli,  they  can  reach  extremely  large  shear  strains  making 
the  nonlinear  response  of  the  tissues  extremely  important  to  characterize.  These  tissues  typically 
are  viscoelastic  in  nature  as  well;  however,  in  this  report  we  do  not  include  rate-dependent 
behavior.  Instead,  we  focus  on  the  fibrous  structures  of  soft  tissues  and  how  this  structural 
anisotropy  affects  the  mechanical  behavior  under  different  loading  mechanisms. 

This  report  introduces  an  anisotropic  constitutive  model  for  modeling  the  intervertebral  discs  of 
the  spine.  Since  this  type  of  model  can  be  used  for  other  soft  biological  tissues,  we  compromise 
between  a  fully  general  model  and  one  that  is  overly  specialized.  Thus,  the  model  is  written  in 
such  a  manner  that  it  can  be  easily  extended  to  capture  the  anisotropy  of  other  biological  tissues, 
such  as  brain  or  skeletal  muscle.  In  this  way,  the  constitutive  model  can  be  thought  of  as  a 
numerical-analytical  tool  for  investigating  the  mechanical  response  of  fibrous  tissue. 

Section  2  introduces  the  background  information  relevant  to  intervertebral  discs.  Section  3 
highlights  the  key  points  in  the  derivation  of  a  transversely  isotropic  hyperelastic  constitutive 
model  with  two  fiber  families.  We  then  verify  the  implementation  of  the  model  in  section  4. 
Section  5  describes  an  algorithm  of  how  we  incorporate  the  fiber  directions  for  an  intervertebral 
disc.  Our  future  applications  of  the  constitutive  model  as  it  will  be  applied  to  the  spine,  and 
could  be  applied  to  the  brain  and  skeletal  muscle,  are  discussed  in  section  6.  Finally,  our 
concluding  remarks  are  presented  in  section  7. 
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2.  Structure  and  Biology  of  the  Spine  and  Intervertebral  Discs 


The  vertebral  column,  also  known  as  the  spine,  is  a  bony  structure  comprised  of  vertebrae  and 
intervertebral  discs,  stacked  alternatively  on  top  of  each  other.  As  seen  in  figure  la,  five  different 
regions  make  up  the  spine:  the  cervical  spine,  the  thoracic  spine,  the  lumbar  spine,  the  sacrum  and 
the  coccyx.  Each  individual  vertebra  is  named  by  referring  to  the  first  letter  of  their  region 
(cervical,  thoracic  or  lumbar),  and,  starting  with  the  most  superior  (highest)  vertebra  in  that 
region,  numbered  consecutively  until  the  most  inferior  (lowest)  vertebra  in  the  region  has  been 
named. 

Between  each  vertebra  is  a  soft  tissue  called  the  intervertebral  disc.  The  intervertebral  discs  play 
a  major  role  in  the  motion  of  the  spine  by  supporting  compressive  forces  experienced  during 
flexion  (bending  forward)  and  extension  (bending  backward),  and  resisting  rotation,  tension,  and 
shear  forces  {11).  An  illustration  of  an  intervertebral  disc  is  shown  in  figure  Ic.  The 
intervertebral  discs  are  made  up  of  two  main  components:  an  inner  gelatinous  region,  known  as 
the  nucleus  pulposus  (grey  region),  and  an  outer  ring  of  fibrosus  cartilage,  known  as  the  annulus 
fibrosus  (white  structures  surrounding  the  grey  region).  The  annulus  fibrosus  is  composed  of 
15-25  concentric  rings  called  lamellae  (four  rings  are  depicted  in  black  in  the  figure)  {12,  13). 
These  rings  consist  of  collagen  fibers  embedded  within  an  extracellular  matrix  (alternating  sets  of 
diagonal  lines).  The  orientation  of  the  fibers  varies  between  adjacent  lamellae,  alternating 
approximately  ±30°  to  the  transverse  plane  of  the  disc. 

Since  the  matrix  of  the  annulus  fibrosus  is  relatively  soft,  the  fibers  are  believed  to  play  a 
prominent  role  in  the  intervertebral  disc’s  response  to  tensile  {14)  and  shear  loading.  At  the 
boundary  between  the  vertebra  and  the  intervertebral  disc  is  a  thin  layer  of  semiporous  bone, 
known  as  the  vertebral  endplate.  The  endplates  of  a  healthy  disc  help  absorb  some  of  the 
pressure  that  results  from  mechanical  loading  of  the  spine  and  prevent  the  nucleus  pulposus  from 
bulging  into  the  adjacent  vertebra  {15). 

Over  the  past  60  years,  there  has  been  a  substantial  effort  to  model  the  spine  and  its  individual 
components  {16).  Computational  models  of  the  spine  provide  researchers  with  an  opportunity  to 
gain  a  more  detailed  understanding  of  the  deformation  stress  state  and  failure  of  the  vertebrae  and 
intervertebral  discs.  Typically,  the  annulus  fibrosus  is  modeled  in  one  of  two  ways:  as  a 
homogeneous  composite  of  the  matrix  and  the  fibers  or  as  an  inhomogeneous  composite  of  the 
matrix  reinforced  by  collagen  fibers  {17).  Shirazi-Adl  found  that  representing  the  annulus 
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fibrosus  as  an  inhomogeneous  composite  gave  more  realistic  results  and  helped  provide  a  better 
understanding  of  the  biomechanical  response  of  the  intervertebral  disc  {14).  Spring  elements, 
truss  elements,  and  rebar  elements  oriented  at  ±30°  to  the  transverse  plane  of  the  intervertebral 
disc  have  all  been  used  by  researchers  to  model  the  fibers  of  the  annulus  fibrosus  (i,  18-20). 


a) 
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ANNULUS  FIBROSUS 


Figure  f.  Bony  anatomy  of  the  spinal  column  (panel  a)  and  a  typical  vertebra  (panel  b). 

Vertebra  are  color-coded  according  to  their  location  classification.  Panel  c  is  an 
illustration  (not  drawn  to  scale)  of  an  intervertebral  disc  showing  the  lamellar 
architecture  of  the  annulus  fibrosus  (white)  which  surrounds  the  nucleus  pulposus 
(grey).  Some  layers  have  been  cut  away  from  the  annulus  fibrosus  to  show  the  fiber 
network  within  the  lamellae.  Note  that  only  four  layers  of  lamellae  are  depicted  in 
the  figure  but  the  annulus  fibrosus  usually  has  15-25  layers.  Collagen  fibers 
(diagonal  lines)  are  oriented  at  ±30°  to  the  transverse  plane  of  the  disc,  with  the 
direction  alternating  between  adjacent  lamellae  in  the  annulus  fibrosus. 
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3.  Transversely  Isotropic  Hyperelastic  Constitutive  Model  With  Two  Fiber 
Families 


Soft  tissues  that  are  comprised  of  fibrous  structures,  such  as  muscles,  ligaments,  tendons, 
intervertebral  discs  and  the  brain,  often  exhibit  strong  anisotropy  along  these  fiber  directions  (27). 
In  this  section,  we  provide  the  assumptions  and  physical  arguments  necessary  for  deriving  a 
constitutive  model  for  representing  a  fibrous  material  as  a  nearly  incompressible,  transversely 
isotropic  hyperelastic  material. 

We  largely  follow  the  work  set  out  by  Weiss  et  al.  (22),  Holzapfel  (23),  Pinsky  et  al.  (24),  and 
Nguyen  and  Boyce  (25).  Instead  of  presenting  a  fully  general  model  and  then  specializing  it  to 
our  application,  we  introduce  simplifying  assumptions  to  tailor  the  derivation  to  our  specific 
application.  We  make  assumptions  appropriate  for  a  nearly  incompressible,  transversely 
isotropic  hyperelastic  material  with  up  to  two  fiber  families  that  do  not  interact  with  one  another, 
nor  the  surrounding  matrix,  and  whose  response  depends  only  on  their  stretch. 

Let  F  be  the  deformation  gradient  describing  the  deformation  of  a  material  relative  to  some 
reference  configuration.  The  polar  decomposition  of  F'  is  given  by 

F  =  RU  =  VR,  (1) 

and  because  Ris  a  properly  orthogonal  rotation  matrix,  the  eigenvalues  of  F  are  the  same  as 
those  for  U  and  V,  the  right  and  left  stretch  tensors.  These  eigenvalues  are  also  the  principal 
stretches,  which  we  denote  as  A*.  We  define  the  right  and  left  Cauchy-Green  deformation 
tensors,  C  and  B,  respectively 


C  =  F^F 

(2a) 

B  =  FF^ 

(2b) 

which  have  the  eigenvalues  \  j. 

The  ratio  of  the  current  specific  volume  to  the  reference  specific  volume  is  the  Jacobian  and  is 
given  by  the  determinant  of  the  deformation  gradient: 


J  =  det  F  . 


(3) 
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The  Jacobian  allows  us  to  define  the  distortional  part  of  the  deformation  gradient: 

F  =  J~^F  (4) 

The  distortional  part  of  the  deformation  gradient  essentially  normalizes  the  volume  changes 
associated  with  the  deformation  and  is  denoted  by  a  bar.  This  can  be  seen  by  taking  the 
determinant  of  F: 


detF  =  det(J-^F)  =  {J-^)det{F)  =  =  1.  (5) 

We  note  that  the  eigenvalues  of  F  are  J~^Xi  =  Aj. 

Analogous  to  the  traditional  approach  in  defining  the  right  and  left  Cauchy-Green  deformation 
tensors,  we  define  the  so-called  modified  right  and  left  Cauchy-Green  deformation  tensors: 

C  =  (6a) 

B  =  FF^.  (6b) 


The  eigenvalues  of  C  and  .B  are  A^. 

We  next  consider  the  invariants  of  the  deformation  tensors  as  they  will  be  important  for  deriving 
our  hyperelastic  constitutive  response  from  an  energy  density  function.  Let  Ji,  I2,  and  J3  denote 
the  first  three  invariants  of  C  and  B: 


Ii  =  trC  =  trB  (7  a) 

J2  =  ^  [iixCf  +  ixC^]  =  i  [iixBf  +  ixB^]  (7b) 

Is  =  det  C  =  det  B  =  (7c) 

and  Ii,  1 2,  and  the  invariants  of  C  and  B: 

1 1  =trC  =  trB  (8a) 

72  =  ^  (trC)VtrC^  (trB)VtrB^  (8b) 

J3  =  det  C  =  det  B  =  1 .  (8c) 


To  incorporate  anisotropy  into  our  description,  we  define  two  fiber  family  directions  in  the 
reference  configuration  ao  and  go,  with  the  property  that  |ao|  =  1  and  |gio|  =  1-  The 


5 


corresponding  deformed  fibers  are  given  by  applying  the  deformation  gradient  to  the  fiber 
direetion  in  the  referenee  eonfiguration  so  that 

a  =  Fao  ,  a  =  Fao  (9a) 

g  =  Fgo  ,  g  =  Fgo  .  (9b) 

The  lengths  of  the  deformed  fiber  families  are 

=  ^{FaoYFao  =  sJao^F^Fao  =  ^ao^Cao  (10) 

and  sinee  |ao|  =  1,  the  fiber  streteh  is 

Aa  =  =  ^Jao^Cao.  (11) 

The  same  arguments  ean  be  made  for  the  seeond  fiber  family,  resulting  in 

Xg  =  V9o^Cgo.  (12) 

Physieal  arguments,  see  for  example,  Weiss  et  al.  (22)  or  Holzapfel  (23),  lead  to  the  eonelusion 
that  the  energy  density  funetion  must  depend  on  an  even  funetion  of  the  fiber  direetions.  Thus, 
one  ean  eonelude  that  the  energy  density  funetion  must  depend  on  the  strueture  tensor  ®  ao 
and  go®  go-  For  notational  simplieity,  let: 

Ao  =  ao®ao,  and  Go  =  go®  go  (13) 

and  in  the  deformed  eonfiguration, 

A  =  a®  a,  and  G  =  g®g.  (14) 

These  additional  tensors  introduee  additional  pseudo-invariants  of  C,  Aq,  and  Gq,  whieh  are 


6 


given  by 


ao^Cao  =  \l 

(15a) 

ao^C^ao 

(15b) 

go^Cgo  =  XI 

(15c) 

go^c^go 

(15d) 

aJCgo 

(15e) 

go  ■ 

(15f) 

Similar 

reasons 


expressions  can  be  derived  for  the  pseudo-invariants  of  the  distortional  tensors,  but  for 


- 


- T-, 


and  Iq, 

(16a) 

(16b) 

The  energy  density  function  0  for  a  hyperelastic  material  is  often  written  in  terms  of  the  right 
Cauchy-Green  deformation  tensor  from  which  the  second  Piola-Kirchhoff  stress  tensor  S  can  be 
determined, 


(17) 


The  energy  density  function  can  equivalently  be  thought  of  as  some  function  of  the  first  three 
invariants  of  either  C  or  B.  Since  the  eigenvalues  of  C  and  C  or  equivalently,  B  and  B,  are 
related  by  J,  the  energy  density  function  can  also  be  expressed  in  terms  of  the  first  three 
invariants  of  either  C  or  B.  Thus,  the  motivation  for  introducing  the  modified  deformation 
tensors  is  that  the  energy  density  can  be  written  as  some  function  of  the  Jacobian  J  and  the 
modified  right  Cauchy-Green  deformation  tensor  C,  i.e., 

<P  =  ^{J,C).  (18) 


Equivalently,  the  energy  density  must  be  some  properly  invariant  function  of  the  nine  invariants 
previously  discussed: 

0  =  Ii,  I2,  Ir-,  I5,  Iq,  I71  Is,  I9) ,  (19) 

where  we’ve  intentionally  isolated  out  the  dependence  of  I3  =  J.  Since  little  is  known  about  the 
actual  constitutive  response  of  these  types  of  materials,  a  common  simplification  introduced  is  to 
assume  that  the  fibers  do  not  interact  with  one  another.  Thus,  we  assume  that  the  mechanical 
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response  is  proportional  to  the  first  three  isotropie  invariants  and  only  the  stretehes  of  the  fiber 
families  (J4,  Jg),  i.e.,  assuming  the  energy  density  is  a  funetion  of  fewer  parameters: 


0  =  0'(J,/i,/2,/4,/6). 


(20) 


An  additional  simplifieation  ean  be  made  when  eonsidering  nearly  ineompressible  materials. 
Typieally  for  these  soft  materials,  the  energy  density  is  deeoupled  into  a  spherieal  part  (relating  to 
pressures  resulting  from  volume  ehange)  and  a  deviatorie  response  (shear  response  independent 
of  volume  ehanges).  This  is  only  approximately  true  for  a  real  material  sinee  there  will  be 
eoupling  of  pressure  and  shear  terms  at  large  deformations.  This  assumption  deeouples  the 
energy  density  as  follows: 

0  =  0s('2)  +  0dev(-^l; -^2, -^4; -^e)  •  (21) 

where  it  ean  be  shown  that  the  pressure  p  is 


P  = 


d4>s 
dJ  ’ 


(22) 


and  the  deviatorie  part  of  the  Cauehy  stress  tensor  dev  T  is 


devT 


(23) 


The  partial  derivative  of  the  energy  density  with  respeet  to  the  distortional  part  of  the  right 
Cauehy-Green  tensor  C  can  be  expanded  using  the  ehain  rule.  This  proeedure  requires  the 
following  additional  partial  derivatives: 


dh 

dC 

(9/2 

W 

(9/4 

w 

dC 


I 


hl-C 


Aq 
Gq  . 


Thus,  the  deviatorie  part  of  the  Cauehy  stress  ean  be  written: 


dev  T  =  —dev 

fj 


—  /  /  (90  ,  ^  (90 


(90  —  ,  (90 


(90 


^  I  -  -^C  +  -^Ao  + 

dh  dh)  dh  dh  dh 


-^T 


(24a) 

(24b) 

(24e) 

(24d) 

(25) 
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—  — T 

Multiplying  through  by  F  on  the  left  and  its  transpose  F  on  the  right  allows  us  to  calculate  the 
deviatoric  part  of  the  Cauchy  stress  in  the  reference  configuration 


dev  T 


\dh  dh)  dh  dh  dh 


(26) 


Following  what  has  been  done  previously  in  the  literature,  we  further  assume  that  the  isotropic 
response  is  that  of  the  well  known  Mooney-Rivlin  model,  so  that  d(f)/dli  =  Cio  and 
(90/(912  =  Coi  are  constants.  We  also  assume  that  the  stress  response  of  both  of  the  fiber  families 
follow  the  same  functional  form,  i.e.,  d(f)/dlQ  =  (90/(914  =  F{I),  so  that  the  deviatoric  part  of  the 
Cauchy  stress  is 


devT 


(C'lo  +  -^iCoi)  B 


Co,B"  +  Fih)A  +  Hh)G 


(27) 


As  is  typical  for  nearly  incompressible  materials,  we  assume  the  spherical  part  of  the  Cauchy 
stress  to  be 

p  =  —Kin  J ,  (28) 

where  k  is  the  bulk  modulus.  Thus,  the  total  Cauchy  stress  is  given  by: 


T  =  Kin  JI  +  —dev 

fj 


(Cio  +  JiCoi)  B  -  CoiB^  +  F{h)A  +  F{h)G 


(29) 


We  have  intentionally  avoided  discussing  the  functional  form  of  the  fiber  response  and  avoided 
writing  down  the  hyperelastic  energy  density  function  explicitly.  The  choice  of  the  fiber  response 
depends  on  the  biological  tissue  being  modeled.  In  the  literature,  collagen  fibers  (25),  as  well  as 
other  fibers,  have  been  modeled  using  an  exponential  function  (22,  24).  This  gives  a  particular 
definition  of  the  fiber  response  F: 


F{I)  =  G  (exp  [A(/  -  1)]  -  1) 


(30) 


where  the  values  of  Ci  and  /3i  can  depend  on  the  fiber  family.  The  fiber  response  function  takes 
the  barred  invariant  /  of  a  fiber  family  and  returns  the  stress  that  results,  thus,  in  practice,  /  in 
equation  30  will  be  either  J4  or  Jg.  While  an  energy  density  function  can  be  written  down  for  the 
case  of  equation  30,  it  may  not  be  possible  for  all  cases  and  functional  forms  of  F.  An  example 
of  this  might  be  including  damage  or  dissipation  to  the  fiber  response.  Section  6.3  discusses  how 
this  model  can  be  extended  to  include  a  prestress  and  section  6.4  discusses  an  active  contractile 
component  that  responds  to  applied  strains. 
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4.  Verification  of  Numerical  Model 


This  section  briefly  covers  various  test  cases  to  illustrate  that  our  constitutive  model  has  been 
implemented  properly  and  behaves  as  expected.  We  conducted  four  simulations  on  a  single 
element  where  we  control  the  input  fiber  direction  and  the  imposed  deformation.  The  first  three 
cases  involve  a  single  fiber  family,  and  the  fourth  incorporates  both  fiber  families.  To  ensure  a  set 
of  rigorous  tests,  we  considered  compression,  extension,  and  shear  cases  for  a  number  of  fiber 
family  orientations.  While  closed-form  solutions  to  each  deformation  were  worked  out  using 
MuPad,  they  are  too  lengthy  to  be  of  any  real  analytical  use.  Instead,  we  compared  our 
simulation  results  directly  against  the  theoretically  predicted  responses  in  various  figures  where 
the  angle  of  the  fiber  families  vary  between  0  and  tt. 

In  the  first  three  tests,  the  initial  fiber  directions  were  taken  so  that  the  fiber  had  no  component  in 
the  X,  y,  or  2:  axis  (corresponding  to  the  planes  YZ,  XZ,  and  XY),  respectively,  i.e., 

ao  =  (0,  cos  6^,  sin  0) ,  ao  =  (cos6',  0,  sin6') ,  or  ao  =  (cos^,  sin^,  0)  (31) 

By  sweeping  6  from  0  to  vr,  we  tested  cases  for  which  the  fiber  direction  was  not  along  a  principal 
axis  of  strain. 

Table  1  lists  the  material  parameters  used  in  this  verification.  In  the  single  fiber  family  cases, 

Cg  =  0.  These  parameters  were  chosen  because  of  their  relevance  to  the  soft  intervertebral  discs 
and  to  illustrate  an  important  issue  regarding  the  sensitivity  of  the  fibers  to  numerical  error  (see 
shear  test  for  a  single  fiber  family).  Additional  tests  (not  shown  here)  further  verified  this  model 
for  additional  choices  of  material  parameters. 


Table  1 .  Summary  of  material  parameters  used  for  the  purpose  of 
verifying  the  numerical  model. 


K  (MPa) 

C7io  (kPa) 

Coi  (kPa) 

Ca(MPa)  /?, 

C73(MPa)  ^g 

7.5 

300 

75 

80  1 

40*  1 

*Cg  =  0  in  the  single  fiber  family  cases 
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4.1  Stretch  and  Compression  Test  for  a  Single  Fiber  Family 


In  these  first  two  tests,  the  element  is  either  stretched  or  compressed  in  the  Z  direction.  The 
results  of  these  test  are  shown  in  figures  2  and  3,  respectively.  In  both  cases,  the  physical 
components  of  the  deformation  gradient  are  given  by 


\F]  = 


1  0  0 

0  1  0 

0  0  q; 


so  that  in  compression  a  <  1,  and  in  tension  a  >  1.  In  both 
Jacobian.  Therefore,  these  tests  also  verify  that  the  pressure 


(32) 


cases,  a  is  also  the  value  of  the 
response  is  implemented  correctly. 


It  is  important  to  note  the  fiber  stretches  that  are  predicted  in 
fourth  pseudo-invariant  is 


u 


al  +  al  + 

2  5 

as 


each  case. 


From  equation  15a,  the 
(33) 


so  that  in  both  the  stretch  and  compression  test 


-  cos  6*^  +  sin  6*^  -  cos  6^^  +  sin  6*^  -  1 

I A  = - ^ - ,  I A  = - ^ - ,  or  I4  =  —  .  (34) 

as  as  as 

This  implies  that  even  for  the  case  where  the  fiber  has  no  2:  component,  the  fiber  response  will  be 
a  constant  and  independent  of  9  since  equation  30  depends  only  on  I4*  Figure  2  shows  the 
results  for  the  single  fiber  family  stretch  tests  where  a  =  1.5.  The  simulation  results  are  shown 
as  the  symbols  (x’s)  and  the  theoretical  response  as  solid  lines.  Note  that  each  symbol  represents 
a  separate  single-element  simulation  with  a  unique  fiber  family  orientation,  so  that  36  simulations 
are  represented  in  each  panel.  Each  component  of  the  Cauchy  stress  T  is  represented  by  a  unique 
color.  The  individual  panels  plot  the  stress  as  it  depends  on  the  angle  9,  which  specifies  the 
undeformed  fiber  directions  given  by  equation  31.  The  setup  for  figure  3  is  the  same,  with  the 
exception  that  a  =  0.5.  Even  at  these  large  deformations,  both  figures  2  and  3  show  excellent 
agreement  between  the  numerical  implementation  of  the  constitutive  model  and  the  theoretical 
predictions  regardless  of  the  fiber  plane  or  angle. 


*  Depending  on  the  available  experimental  data  and  or  the  qualitative  features  desired  from  the  model,  equation  30 
could  be  rewritten  to  depend  on  I4  instead  of  I4  with  minimal  alterations  to  the  code. 
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Figure  2.  Single  fiber  family  stretch  test.  Simulation  (symbols)  comparison  against  theoretical  (solid 

lines)  for  the  components  of  the  Cauchy  stress  Txx  (red),  Tyy  (blue),  T^z  (black),  Txy  (cyan), 
Tyz  (magenta),  and  T^x  (green)  for  three  sets  of  fiber  orientation  vectors  oq.  Panel  a  shows  the 
stress  response  of  a  single  fiber  family  with  initial  direction  vector  in  the  Y Z-plane.  Similarly, 
panels  b  and  c  show  a  fiber  family  with  initial  direction  vector  in  the  XZ-plane,  and  XY -plane, 
respectively. 


Figure  3.  Single  fiber  family  compression  test.  Simulation  (symbols)  comparison  against  theoretical 
(solid  lines)  for  the  components  of  the  Cauchy  stress  Txx  (red),  Tyy  (blue),  T^z  (black),  Txy 
(cyan),  Tyz  (magenta),  and  Tzx  (green)  for  three  sets  of  fiber  orientation  vectors  a^.  Panel  a 
shows  the  stress  response  of  a  single  fiber  family  with  initial  direction  vector  in  the  Y Z-plane. 
Similarly,  panels  b  and  c  show  a  fiber  family  with  initial  direction  vector  in  the  2fZ-plane,  and 
XY -plane,  respectively. 
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4.2  Shear  Test  for  a  Single  Fiber  Family 

The  third  case  considered  for  the  single  fiber  family  was  a  shearing  in  the  V  direction  in  the 
Z-plane  (figure  4).  In  this  case,  the  physical  components  of  the  deformation  gradient  are  given  by 

■  1  0  0  ■ 

[F]=  Ola  .  (35) 

_  0  0  1  _ 

Thus  the  fourth  pseudo-invariant  takes  on  the  form 

1 4.  =  1  T  2q;  (Xy  (Xz  -f  Oi  Qj^  (36) 

so  that  in  each  of  the  cases  in  equation  3 1 : 

J4  =  1  +  2q;  cos  6*  sin  0  +  sin  ,  1^  =  1  +  a^  sin  ,  or  14  =  1 .  (37) 

Figure  4  compares  the  simulation  results  for  the  single  fiber  family  in  shear.  Panels  a  and  b  show 
excellent  agreement  between  the  theoretically  predicted  values  and  the  simulation  results. 

Panel  c,  however,  shows  some  noteworthy  deviations  from  the  predicted  theoretical  curves.  The 
deviation  from  the  theoretical  curves  is  an  important  issue  to  highlight  and  is  entirely  due  to  the 
small  numerical  error  that  SIERRA  introduces  in  the  rotation  matrices  during  integration  steps. 
Note  that  the  scale  bars  in  panel  c  are  in  kPa,  while  the  other  two  panels  are  in  MPa.  Also  note 
that  for  this  specific  example  (see  table  1),  the  choice  of  moduli  place  the  isotropic  response  three 
orders  in  magnitude  smaller  than  the  fiber  response.  This  particular  test  is  an  extreme  case  where 
the  deformation  leaves  the  fibers  unchanged  so  that  1 4  should  be  unity  and  that  by  equation  30, 
the  fiber  response  should  be  0,  leaving  only  the  isotropic  stress  response.  Upon  closer  inspection, 
the  rotation  matrices  calculated  by  SIERRA  had  small  errors  when  compared  to  the  theoretically 
predicted  values.  In  the  pure  shear  case  described  previously,  components  of  the  left-Stretch 
tensor  V  and  the  rotation  tensor  R  should  be  as  follows: 
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Figure  4.  Single  fiber  family  shear  in  Y-direction  test.  Simulation  (symbols)  comparison  against 

theoretical  (solid  lines)  for  the  components  of  the  Cauchy  stress  T^x  (red),  Tyy  (blue),  Tzz 
(black),  Txy  (cyan),  Ty^  (magenta),  and  T^x  (green)  for  three  sets  of  fiber  orientation  vectors 
ao-  Panel  a  shows  the  stress  response  of  a  single  fiber  family  with  initial  direction  vector  in  the 
Y Z-plane.  Similarly,  panels  b  and  c  show  a  fiber  family  with  initial  direction  vector  in  the 
YZ-plane,  and  XY-plane,  respectively. 


In  practice,  however,  we  found  even  with  striet  demands  on  eonvergenee  eriteria  (in  the  implieit 
ease)  or  small  time  steps  (in  the  explieit  ease),  the  yy,  yz,  zy,  and  2:2;  eomponents  of  the  aetual 
tensors  deviated  from  the  theoretieally  predieted  values 


[^S 


1 

0 

0 


2+a-^ 

\/4+a 

a 

y/  4+a^ 


—  € 


+  6/2 


(39) 


1 

0 

0 


V4+C 


\/4+, 


0 

I  -  e/10 

+  e/2 


0 

V4+a2  “ 


(40) 


We  estimated  e  from  a  few  numerieal  simulations  of  the  shear  tests  and  found  that  its  value 
typieally  was  small,  e  ~  10“^.  The  produet  of  V  and  R  is  F,  the  deformation  gradient,  whieh  to 
first  order  in  e  is 


1  0  0 
0  l-efyy{a)  a-efy^{a) 
0  -efzy{a)  1  +  efzz{a) 


(41) 
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where  fij{a)  are  functions  of  the  deformation.  These  functions  were  close  to  unity  and  their 
functional  form  is  suppressed  since  they  are  not  essential  to  the  arguments  that  follow.  Using  this 
form  of  the  deformation  gradient,  we  can  obtain  an  approximate  expression  for  the  fourth 
pseudo-invariant  to  the  first  order  in  e: 


h-l-ef{a,9).  (42) 

Inserting  this  expression  into  equation  30  gives 

J'(/4)  =  U3  [exp{-/3aef{a,  6))  -  1]  -eC^/Safia,  9) .  (43) 

Using  the  values  from  table  1  and  e  =  10“"^  gives  a  fiber  response,  on  the  order  of  10  kPa,  whose 
magnitude  will  depend  further  on  the  deformation  and  the  angle  of  the  fiber.  This  is  all  of  the 
same  order  as  the  error  in  figure  4c.  Thus,  the  choice  of  an  exponential  function  can  make  the 
simulation  very  sensitive  to  numerical  error. 

4.3  Compression  Test  With  Two  Fiber  Families 

The  previous  simulations  were  performed  for  single  fiber  families.  The  final  verification  of  our 
implementation  is  a  two  fiber  family  test,  in  which  we  verify  that  the  fiber  families  behave 
properly  when  both  families  are  used.  For  this  test,  the  fiber  families  are  represented  by  the 
vectors: 

ao  =  (cos6',0,sin6*) ,  and  gro  =  (0,  sin^,  cos  0) .  (44) 

This  choice  of  fiber  directions  can  be  either  orthogonal  or  non-orthogonal  (ao  ■  go  =  sin  9  cos  9) 
depending  on  9.  Figure  5  compares  the  simulation  versus  theoretical  results  for  a  compression 
test.  By  varying  9,  our  choice  of  ao  and  Qq  sweep  the  fiber  families  through  the  XZ-  and 
Y Z-planes,  respectively.  This  figure  also  shows  excellent  agreement  between  the  theoretically 
predicted  response  and  the  simulation  results. 
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Figure  5.  Two  fiber  families  compression  test.  Simulation  (symbols)  comparison  against  theoretical 

(solid  lines)  for  the  components  of  the  Cauchy  stress  T^x  (red),  Tyy  (blue),  Tzz  (black),  Txy 
(cyan),  Tyz  (magenta),  and  T^x  (green)  for  two  fiber  family  orientation  vectors  oq  and  Qq. 


5.  Determining  the  Fiber  Directions  for  an  Intervertebral  Disc 


This  section  discusses  the  overall  approach  of  how  we  approximate  the  fiber  family  directions 
within  intervertebral  discs.  We  begin  with  a  more  mathematical  description  of  the  fiber 
orientations  within  a  simplified  intervertebral  disc  and  then  discuss  details  of  an  algorithm  to 
handle  some  more  general  geometries. 

Our  constitutive  model  is  designed  to  read  in  a  file  that  lists  the  fiber  family  directions  as  they 
change  in  space.  Although  this  model  was  designed  initially  for  intervertebral  discs,  it  can  in 
theory  be  extended  to  handle  a  number  of  other  materials  that  have  the  same  feature  of  one  or  two 
fiber  directions,  e.g.,  collagen  fibers  of  the  cornea,  striated  muscle  fibers  in  skeletal  muscles, 
multiple  axonal  directions  within  the  brain.  In  each  case,  the  choice  of  the  fiber  direction,  or  the 
manner  in  which  it  is  assigned  to  an  element,  could  vary  drastically.  Since  we  typically  keep  the 
same  mesh  from  one  simulation  to  another,  we  separated  the  calculation  of  the  fiber  directions  per 
element  into  a  preprocessing  step  handled  in  MATLAB.  This  choice  avoided  repeated  upfront 
costs  associated  with  determining  where  fibers  were  with  respect  to  a  mesh,  and  enables  rapid 
changes  to  be  made  without  requiring  a  recompile  of  SIERRA. 
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As  discussed  in  section  2,  the  intervertebral  discs  are  reinforced  by  collagen  fibers  in  the  annulus 
fibrosus.  Figure  Ic  shows  that  these  fibers  typically  have  a  regular  arrangement  within  the 
lamellae  of  the  annulus  fibrosus.  To  describe  these  directions,  we  start  by  approximating  the 
intervertebral  disc  as  a  cylinder,  as  shown  in  figure  6.  Next,  we  consider  a  point  on  the  surface  of 
the  cylinder,  in  the  figure  this  is  given  by  some  vector  r  in  the  reference  frame.  This  surface 
corresponds  to  a  single  lamellae  ring.  A  surface  normal  can  be  defined  on  the  cylinder  as  n  for 
which  there  is  a  tangent  plane.  In  this  example,  we  set  the  tangent  vector  t  to  be  parallel  with  the 
Z-axis,  however,  in  general  this  will  not  be  the  case  for  an  intervertebral  disc  since  it  will  be  taken 
to  coincide  with  the  normal  to  the  transverse  plane  of  the  intervertebral  disc.  By  definition,  the 
binormal  vector  is  b  =  t  x  h.  According  to  the  experimental  literature  (12,  19),  the  two  fiber 
family  orientations  and  are  perpendicular  to  the  surface  normal  vector  n  and  make  an  angle 
9  with  the  tangent  vector,  i.e.. 


and 


tto  ■  n  =  0  ,  ttQ  ■  i  =  cos  9a 
gQ-h  =  0,  g^ -i  =  cos  9  g 


(45) 

(46) 


(47) 


Figure  6.  Local  coordinate  system  for  an  intervertebral  disc. 
An  idealized  intervertebral  disc  in  the  reference 
coordinate  system.  The  point  r  is  on  the  surface  of 
the  cylinder  with  corresponding  normal  vector  n, 
chosen  tangent  vector  t,  and  binormal  vector  b.  The 
vectors  uq  and  Qq  for  this  point  are  also  shown. 
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If,  in  addition  to  the  cylinder  shown  in  figure  6,  a  second  concentric  cylinder  with  a  smaller  radius 
is  added,  it  would  share  the  same  fiber  orientations.  The  experimental  literature  describes  a 
similar  concentric  structure  to  the  orientation  of  the  fiber  families  in  the  lamellae  rings  within  the 
annulus  fibrosus.  A  more  generic  notion  of  this  concept  will  be  used  later  when  we  describe  how 
we  determine  fiber  directions  for  a  specific  annulus  fibrosus  geometry.  Thus,  if  one  can 
determine  a  surface  normal  and  the  tangent  vector  that  is  parallel  to  the  intervertebral  disc’s 
vertical  axis  (normal  to  the  transverse  plane  of  the  intervertebral  disc),  then  one  can  calculate  the 
fiber  family  orientations  by  applying  the  appropriate  rotations. 

To  extend  the  application  of  these  mathematical  concepts  to  a  more  complicated  geometry  using  a 
semi-automated  approach,  the  overall  procedure  is  split  between  two  programs:  a  command  line 
tool  and  MATLAB. 

The  first  step  is  handled  through  a  collection  of  scripts  and  a  tool  provided  in  the  SEACAS 
toolbox  that  is  included  with  SIERRA,  namely  “GROPE”.  This  command  allows  one  to 
manipulate  and  survey  the  mesh  through  the  command  line  and  send  results  to  an  ASCII  file.  We 
used  GROPE  to  extract  the  centroid  location  of  each  element  within  a  given  annulus  fibrosus  and 
saved  it  to  a  delimited  file.  We  use  this  centroid  data  as  a  point-cloud  approximation  for  the 
actual  geometry  that  we  manipulate  in  MATEAB.  A  mesh  of  the  spinal  segment  E3E4E5  is  shown 
in  figure  7a  and  a  sagittal  view  of  the  corresponding  point-cloud  for  a  single  annulus  fibrosus  is 
shown  in  figure  7b. 
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Figure  7.  Determining  fiber  directions  in  the  spine.  Panel  a  is  the  meshed  L3L4L5  segment  of  our  spine 
model.  In  this  drawing,  the  anterior  (forward  facing)  side  corresponds  to  the  F-axis  and  the 
vertical  is  Z.  The  centroids  of  all  the  elements  of  the  beige  intervertebral  disc  are  shown  as 
blue  dots  in  panel  b.  A  linear  fit  to  the  top  layer  of  centroids  (cyan)  gives  the  slope  of  that 
layer.  The  angle  normal  to  the  plane  of  the  intervertebral  disc  can  be  determined  from  this 
slope  (red  arrow)  relative  to  the  Z-axis  (black  arrow).  This  top  layer  is  shown  after  it  is  rotated 
and  projected  into  the  XY -plane  in  panel  c.  A  convex  hull  algorithm  determines  the 
outermost  ring  of  centroids  (cyan  open  circles)  from  which  a  parameterization  yields  the 
binormal  vector  h  (red  arrows).  Using  the  angle  normal  to  the  plane  of  the  intervertebral  disc 
as  the  surface  tangent  vector  f  and  the  binormal  vector  from  the  convex  hull  h,  the  fiber 
orientations  can  be  determined.  Panel  d  shows  a  close-up  of  a  wireframe  of  the  original  mesh 
where  the  calculated  two  fiber  families  are  shown  in  red  and  green. 


The  algorithm  developed  in  MATLAB  is  broken  into  several  steps.  It  is  applied  to  a  set  of 
intervertebral  discs,  but  we  only  discuss  its  application  to  a  single  intervertebral  disc. 

The  first  step  of  the  MATLAB  code  was  to  simplify  the  three-dimensional  point  cloud  down  to  a 
stack  of  two-dimensional  points.  This  step  is  guided  by  user  input  so  that  the  centroids  of 
elements  that  all  belonged  to  the  same  layer  of  the  mesh  could  be  selected.  This  simplifying  step 
was  only  possible  because  a  clear  sweep  direction  of  the  mesh  could  be  defined.  Figure  7b  shows 
that  a  view  of  the  XZ-plane  provides  enough  space  between  layers  to  differentiate  them.  Each 
layer  of  the  intervertebral  disc  is  a  collection  of  three-dimensional  points  that  approximately 
resided  in  a  plane.  At  this  step  we  made  use  of  another  symmetry  of  the  intervertebral  discs,  the 
L-R  symmetry,  which  in  this  case  is  oriented  with  the  Y -axis.  A  linear  fit  of  the  {X,  Z) 
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coordinate  points  within  a  single  layer  gave  an  approximate  slope  of  that  plane  (see  figure  7b). 
This  slope  is  then  used  to  obtain  a  normal  vector  (red  arrow  in  panel  b),  which  typically  fell 
within  20°  of  the  Z-axis  (black  arrow).  It  is  important  to  note  that  the  term  “normal  vector” 
refers  to  the  fact  that  the  vector  is  normal  to  the  plane  of  the  single  layer  of  element  centroids. 
However,  in  our  formulation,  this  direction  actually  represents  the  tangent  vector  t  from  figure  6 
since  it  is  tangent  to  the  surface  of  the  annulus  fibrosus.  Using  the  slope,  the  layer  is  rotated  so 
that  the  coordinates  of  the  centroids  primarily  fell  in  the  XY -plane.  With  the  tangent  vector  i 
approximately  known  for  this  layer,  it  remains  to  determine  the  outward  normal  n  or  the  binormal 
vector  b.  We  note  that  our  procedure  assumes  that  the  average  slope  of  a  layer  of  elements 
provides  a  good  estimate  for  determining  the  tangent  vector  t.  However,  this  will  not  be  the  case 
if  the  intervertebral  disc  exhibits  large  barreling.  In  this  extreme  case,  additional  considerations 
might  be  necessary  to  approximate  n  or  t. 

Figure  7c  shows  the  result  of  rotating  and  projecting  the  centroids  of  the  top  layer  of  panel  b  to 
the  XY -plane.  Since  the  centroids  of  the  annulus  fibrosus  essentially  form  a  set  of  elliptical  rings 
in  this  plane,  we  use  a  convex  hull  algorithm  to  single  out  the  centroids  that  belong  to  the 
outermost  ring  (cyan  open  circles,  panel  c).  Using  the  coordinates  of  the  outermost  ring,  a 
parametric  description  of  the  elliptical  ring  is  formed,  i.e.,  (X (s),  U(s)).  The  tangent  line  of  the 
parametric  curve  given  by  (X(s),  Y (s))  then  represents  the  binormal  vector  b,  the  result  of  this 
calculation  is  shown  as  the  red  arrows.  This  procedure  can  be  repeated  by  eliminating  the  outer 
ring  of  points  and  reapplying  a  convex  hull  algorithm  to  identify  the  next  layer  of  centroids.  We 
emphasize  that  this  algorithm  heavily  depends  on  the  regularity  of  our  mesh  and  that  other 
algorithms  may  be  necessary  for  more  complicated  meshes  or  geometries. 

The  fiber  orientations  for  the  full  intervertebral  disc  in  the  lab  frame  can  be  found  by  using  a 
specific  fiber  family  angle  (±9  in  figure  6)  given  relative  to  the  tangent  vector  t  (with  associated  b 
and  n).  This  is  done  by  working  in  the  reference  frame  of  a  single  layer  of  elements  in  the 
intervertebral  disc,  i.e.,  the  slope  determined  from  the  single  layer  is  used  to  obtain  X^isc  and  the 
perpendicular  (which  corresponds  to  t)  is  our  Z^isc-  We  then  assign  a  fiber  vector  to  coincide  with 
^disc,  which  in  this  frame  is  (0,0,1).  The  first  step  is  rotating  the  fiber  vector  about  the  Idisc-axis 
so  that  it  makes  an  angle  9  with  the  Z^isc  vector  (or  —9  depending  on  the  fiber  family).  The  next 
rotation  accounts  for  the  measured  binormal  vector  b  for  a  given  element  by  rotating  the  resultant 
vector  about  the  Zdisc-axis  by  the  angles  represented  by  the  red  arrows  in  figure  7c.  The  final 
rotation  is  to  give  the  single  layer  of  elements  the  appropriate  slope  determined  earlier  by  the 
algorithm.  This  is  accomplished  by  rotating  the  resultant  vector  about  the  Idisc-axis  by  the  angle 
between  the  Zdisc-axis  and  the  Z-axis,  an  angle  that  was  typically  less  than  20°.  The  fiber  vectors 
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are  then  translated  to  the  element  centroid  location,  and  then  translated  again  to  return  the  annulus 
fibrosus  to  coincide  with  the  reference  frame. 


In  practice,  this  procedure  works  well  for  objects  whose  cross  section  is  convex.  The  results  for 
two  intervertebral  discs  are  shown  in  figure  7d,  which  shows  a  magnified  wireframe  view  of  the 
original  mesh.  The  fiber  families  are  represented  by  two  sets  of  vectors  ao  (red)  and  Qq  (green). 
While  this  procedure  might  only  approximate  experimental  data,  the  strength  of  our  approach  is 
that  we  can  substitute  this  preprocessing  step  with  actual  experimental  data  in  the  future. 


6.  Future  Applications  of  the  Model 


This  section  describes  some  applications  of  our  transversely  isotropic  hyperelastic  constitutive 
model  with  two  fiber  families.  The  first  example,  discussed  in  section  6. 1  presents  a  preliminary 
result  in  modeling  an  intervertebral  disc.  The  second  example,  section  6.2,  explores  the 
possibility  of  applying  this  model  to  the  brain  as  a  continuation  of  previous  research  which 
considered  a  transversely  isotropic  hyperelastic  constitutive  model  with  a  single  fiber  family.  The 
last  two  examples  explore  how  the  model  could  be  generalized  to  incorporate  a  prestress 
(section  6.3)  and  an  active-contractile  element  similar  to  skeletal  muscle  (section  6.4). 

Meshes  were  generated  in  Cubit  (V13. 1 ;  Sandia  National  Laboratory).  Simulations  were 
performed  using  SIERRA/SolidMechanics  (Adagio/Presto  4.28;  Sandia  National  Laboratory). 
Adagio  is  an  implicit,  nonlinear  preconditioned  conjugate  gradient  solver  and  Presto  is  an  explicit 
solver.  Postprocessing  of  simulation  results  was  carried  out  in  Para  View  (V3.14.0;  Kitware)  and 
MATLAB  (The  Math  Works,  Natick,  MA).  Preprocessing  the  fiber  directions  was  performed 
using  GROPE  (SEACAS  Toolkit,  Sandia  National  Eaboratory),  and  MATEAB.  Additional 
theoretical  calculations  were  performed  using  MuPad  an  application  package  of  MATEAB. 
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6.1  Modeling  Intervertebral  Discs 


This  section  briefly  covers  an  example  application  of  our  constitutive  model  applied  to  the 
intervertebral  discs  and  compares  it  against  an  isotropic  Mooney-Rivlin  material.  We  used  the 
algorithm  described  in  section  5  to  approximate  fiber  family  directions  for  a  single  intervertebral 
disc  surrounded  by  two  vertebrae.  In  this  geometry,  the  intervertebral  disc  is  tied  to  the  vertebrae. 
We  quasi-statically  imposed  a  displacement  to  the  top  vertebra  equating  to  a  1%  compressive 
strain,  and  fixed  the  bottom  vertebra.  This  produced  a  slight  barreling,  or  radial  bulging  of  the 
intervertebral  disc. 

The  results  of  our  simulations  are  shown  in  figure  8.  For  the  transversely  isotropic  hyperelastic 
model  with  two  fiber  families,  we  used  the  material  parameters  from  table  1.  The  same  values 
were  used  for  the  Mooney-Rivlin  material  with  the  exception  that  the  fibers  were  turned  off,  i.e., 
Ca  =  Cg  =  0.  For  both  material  models,  the  annulus  fibrosus  was  assumed  to  have  a  density  of 
1200  kg/m^  and  the  nucleus  pulposus  had  a  density  of  1000  kg/m^. 
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Transversely  Isotropic  Model  Mooney- Rivlin  Model 

With  Two  Fiber  Families 


Figure  8.  Comparison  of  constitutive  models.  Panels  a,  c,  e,  g,  i,  k,  and  m 

summarize  the  simulation  results  for  our  transversely  isotropic  hyperelastic 
constitutive  model  with  two  fiber  families.  Panels  b,  d,  f,  h,  j,  1,  and  n  are 
the  results  of  that  same  material  without  a  fiber  response.  In  both  cases, 
the  specimen  is  subjected  to  a  1%  compressive  strain.  The  resultant 
pressure,  and  six  components  of  stress  are  shown  in  false  color.  The 
orientation  of  the  intervertebral  discs  relative  to  a  lab  coordinate  system  is 
shown  in  panel  n  and  the  color  bar  applies  to  all  panels. 


23 


Figure  8  panels  a,  c,  e,  g,  i,  k,  and  m  show  the  pressure  and  six  components  of  the  Cauchy  stress 
T  for  the  intervertebral  disc  with  the  transversely  isotropic  hyperelastic  constitutive  model  with 
two  fiber  families.  The  panels  b,  d,  f,  h,  j,  1  and  n  show  the  corresponding  pressure  and  Cauchy 
stress  for  the  Mooney-Rivlin  intervertebral  disc.  The  false  color  of  each  of  the  intervertebral 
discs  is  set  to  the  same  color  scale  so  that  one  can  directly  compare  the  results  of  the  two  models. 

As  expected,  the  largest  stresses  were  seen  in  the  annulus  fibrosus  in  the  transversely  isotropic 
hyperelastic  constitutive  model  with  two  fiber  families,  particularly  in  the  normal  stresses  T^x  and 
Tyy  (figure  8c  and  e),  but  also  in  the  shear  stresses,  particularly  Txy  (figure  8i).  The  stresses  are 
larger  due  to  the  fiber  response  function,  equation  30.  For  our  choice  of  (3,  the  exponential 
function  is  extremely  sensitive  to  any  stretch  of  the  fiber  family.  For  a  cylinder  under  uniaxial 
compression,  where  the  material  is  allowed  to  expand  radially,  we  would  not  expect  shear  stresses 
in  the  Mooney-Rivlin  material.  However,  since  the  intervertebral  discs  were  tied  to  the  vertebrae, 
the  deformation  exhibited  a  slight  barreling,  and  the  original  geometry  became  irregular  (see 
figure  7a),  small  shear  stresses  developed  even  in  the  Mooney-Rivlin  material,  e.g.,  figure  8 
panels  j,  1,  and  n.  These  stresses  were  smaller  than  those  that  develop  in  the  transversely 
isotropic  hyperelastic  constitutive  model  with  two  fiber  families  (compare  panels  i  and  j,  panels  k 
and  1,  or  panels  m  and  n).  The  pressures  are  larger  in  the  transversely  isotropic  hyperelastic 
constitutive  model  with  two  fiber  families,  most  notably  in  the  nucleus  pulposus  (figure  8a). 

These  observations  from  this  preliminary  quasi-static  simulation  appear  to  support  the  hypothesis 
that  during  dynamic  axial  compression  of  the  intervertebral  disc,  the  pressure  inside  the  nucleus 
pulposus  increases,  forcing  the  endplates  to  bulge  into  the  cancellous  core  of  the  adjacent 
vertebra,  potentially  causing  a  burst  fracture  (26). 

6.2  Modeling  Brain  Tissue 

This  section  briefly  outlines  a  potential  application  of  our  constitutive  model  to  brain  tissue. 
Previously,  a  transversely  isotropic  hyperelastic  constitutive  model  with  a  single  fiber  family  was 
developed  by  Kraft  and  Dagro  (27)  and  later  used  by  McKee  et  al.  (28).  In  this  model,  diffusion 
spectrum  imaging  (DSI)  was  used  to  provide  the  constitutive  model  with  the  corresponding  fiber 
directions.  The  experimental  data,  however,  often  had  more  than  one  unique  fiber  direction  for  a 
given  voxel  of  data.  The  authors’  solution  was  to  provide  an  average  over  fiber  directions  as  a 
first  approximation  with  the  intent  to  extend  their  model  in  the  future  to  handle  additional  fiber 
families.  Using  some  raw  voxel  data  provided  by  these  authors  (27,  28),  we  developed  a  simple 
method  to  extract  two  fiber  directions.  Figure  9  presents  some  example  results  of  this  procedure. 
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Figure  9.  An  approach  to  extract  two  fiber  families  from  diffusion  spectrum 
imaging  (DSI).  Panels  a  and  c  are  the  raw  data  from  DSI  from  two 
exemplary  voxels.  The  color  scale  represents  the  relative  strength  of  the 
specific  direcfion  sampled.  The  firsl  is  a  case  where  fhere  is  only  a  single 
material  direcfion  and  fhe  second  is  a  case  where  fhere  are  fwo  material 
direcfions.  Panels  b  and  d  show  fhe  resulf  of  applying  our  mefhod  fo 
approximafe  fhe  direcfions  from  fhe  raw  dafa. 


The  measurement  reeeived  from  DSI  is  a  set  of  voxels  that  have  a  set  of  sampled  direetions  and 
the  assoeiated  image  intensity  in  that  direetion.  These  data  are  often  represented  by  a  set  of  unit 
veetors  rti  whose  eorresponding  spherieal  eoordinates  are  (p,  6,  cj)),  with: 

P  =  ^  ^  ^  cos“^  (uz)  .  (48) 

In  addition  to  these  sample  direetions  is  a  weight  Wi,  i.e.,  the  magnitude  of  the  measurement 
along  that  sampled  direetion. 

Figure  9a  shows  the  data  obtained  from  a  single  voxel  of  DSI.  Eaeh  filled  eirele  in  this  panel 
represents  a  sampled  direetion  rii  (with  some  Oi  and  (pi).  Its  eolor  represents  the  weight  Wi  or 
intensity  of  the  DSI  measurement  in  the  n*  direetion.  This  partieular  voxel  data  shows  a  ease 
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where  a  single  fiber  family  is  present.  Figure  9c  shows  similar  data  for  a  case  where  two  distinct 
fiber  directions  can  be  seen.  These  two  cases  were  chosen  to  illustrate  the  importance  of 
improving  the  average  weighting  to  one  that  can  identify  two  fiber  directions. 


We  followed  the  same  procedure  for  both  cases  to  extract  the  relevant  key  directions  from  the  raw 
data.  First,  we  linearly  interpolated  the  DSI  data  for  a  single  voxel  and  ran  a  nonlinear  search 
algorithm  (“fminsearch”  in  MATLAB)  to  find  the  maximum  intensity  and  its  corresponding  pair 
of  angles  {9o,  (po)  .  After  finding  the  largest  peak,  we  used  the  Lorentzian  function: 


a 


a  [{9  -  6*0)2  +  (0  -  (poY]  +  1  ’ 


(49) 


centered  at  the  determined  values  6*0  and  0o  to  subtract  off  this  peak.  For  this  part  of  the 
algorithm,  the  DSI  data  are  assumed  to  be  tt  periodic  in  both  9  and  0.  We  chose  a  from  the  peak 
value  and  used  an  empirically  determined  constant  value  of  a  =  0.3.  After  performing  this 
subtraction,  any  negative  weights  were  set  to  zero,  so  that  upon  repeating  the  entire  procedure,  the 
second  largest  peak  could  be  determined  and  removed.  Determining  the  two  largest  peaks  in  the 
9  and  0  interpolated  space  is  equivalent  to  finding  the  two  sets  of  angles  that  best  correspond  to 
the  two  fiber  family  vectors  and  g^. 

Applying  this  procedure  to  the  DSI  data  in  figure  9a  resulted  in  a  calculated  peak  at  {9a=l-51, 
0a=1.57).  The  Lorentzian  fit  for  this  first  case  is  shown  in  figure  9b.  Applying  the  same 
procedure  to  the  DSI  data  in  figure  9c  resulted  in  two  calculated  peaks,  (6*a=3.03,  0a=1.57)  and 
(6*3=1.30,  03=1.57).  The  sum  of  the  two  Lorentzian  functions  for  the  second  case  is  shown  in 
figure  9d.  The  sets  of  angles  given  by  6*  and  0  can  then  be  used  to  determine  normal  vectors  for 
the  two  fiber  families  and  g^.  It  is  important  to  note  that  the  approach  used  in  Kraft  and 
Dagro  (27)  was  to  take  the  weighted  average,  i.e.. 


ao  = 


y^Xwi{nx)i , Wi{ny)i ,  Wi{n^)i ) , 


(50) 


and  enforce  that  |ao|  =  1.  To  illustrate  a  potential  pitfall  of  using  a  weighted  average,  we  also 
calculated  the  corresponding  weighted  average  angles  (6*a,0a)weighted  for  these  two  cases.  Using 
the  method  from  Kraft  and  Dagro  (27),  we  calculated  (6*^,  0a)weighted  =  (3.03, 1.57)  for  the  first 
case  and  (6*^,  0a)weighted  =  (0.04, 1.56)  for  the  second  case.  In  the  first  case,  the  weighted  average 
does  not  return  a  valid  representation  of  the  DSI  data,  since  the  average  is  overwhelmed  by  the 
low-intensity  data.  Restricting  the  average  to  consider  only  the  stronger  weights  (Wi  >  0.7) 
resulted  in  {9a,  0a)weighted  =  (2.15, 1.57),  which  is  a  much  more  reasonable  result.  The  potential 
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issue  with  the  method  from  Kraft  and  Dagro  (27)  is  that  the  data  set  is  represented  by  its  mean 
instead  of  its  most  probable  value.  Additionally,  in  the  cases  where  two  fiber  families  are  present, 
the  weighted  average  cannot  pick  up  on  the  second  peak  that  is  clearly  visible  in  the  raw  data. 

One  can  attempt  to  avoid  the  case  of  two  distinct  fiber  families  by  refining  the  mesh.  This  was 
likely  the  approach  in  Kraft  and  Dagro  (27),  where  a  voxel  of  DSI  data  was  approximately  the 
same  size  as  an  element.  However,  if  a  coarser  mesh  size  is  used,  then  the  number  of  DSI  voxels 
contained  within  an  element  would  increase.  The  odds  of  finding  multiple  fiber  families  within 
an  element  would  then  also  increase.  The  use  of  a  coarse  mesh  and  two  fiber  families  enables 
one  to  capture  the  mechanical  features  that  the  anisotropy  introduces  and  also  benefits  from  a 
larger  time  step  in  an  explicit  code. 

6.3  Incorporating  Fiber  Prestresses 

This  section  describes  how  one  could  accomplish  adding  a  prestress  to  the  fiber  families  through 
their  stretches.  Solvers,  like  SIERRA,  typically  assume  that  the  configuration  read  into  the 
program  (in  the  form  of  the  mesh)  is  the  undeformed  reference  configuration.  In  the  human 
body,  however,  the  geometry  of  an  intervertebral  disc  or  a  muscle  captured  from  imaging,  may 
not  be  the  “elastically  neutral”  or  stress-free  reference  configuration.  It  is  reasonable  to  assume 
that  there  would  be  stresses  present  at  the  very  beginning  of  the  simulation,  e.g.,  intervertebral 
discs  supporting  the  weight  of  the  upper  body,  or  tensed  or  stretched  muscles. 

As  before,  let  F  be  the  deformation  gradient  from  a  reference  configuration  which  has  no  stress. 
Let  lower  case  coordinates  denote  the  current  configuration  (x,  y,  z),  and  upper  case  coordinates 
(X,  F,  Z)  denote  the  reference.  Then  the  deformation  gradient  is  a  linear  transformation 
F  :  {X,  Y,  Z)  — )■  {x,  y,  z).  Assume  also  that  there  is  a  second  reference  configuration  where  the 
stresses  are  not  necessarily  zero.  The  deformations  present  in  this  configuration  would 
correspond  to  the  current  mesh  that  is  sent  to  a  simulation.  Let  the  coordinates  of  this 
configuration  be  represented  by  uppercase  with  primes  (X'Y'Z').  We  can  think  of  a  deformation 
gradient  that  takes  us  from  the  reference  configuration  to  the  configuration  read  into  a  simulation 
as,  F'  :  {X,  Y,  Z)  — )■  (X',  F',  Z').  The  simulation  calculates  a  rotation  matrix  and  left  stretch 
tensor  relative  to  the  primed  configuration  {X'Y' Z'),  and  so  we  can  also  think  of  a  third 
deformation  gradient  Fg  ;  (X',  F',  Z')  — )■  {x,  y,  z)  (subscript  s  for  simulation).  These 
deformation  gradients  are  related  by  the  following  expression: 

F  =  FgF'.  (51) 

Now  we  consider  a  fiber  family  whose  direction  in  the  reference  configuration  is  given  by  ao. 
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where  |ao|  =  1.  We  can  again  think  of  the  deformed  fiber  in  the  current  configuration: 

a  =  Fao  =  FgF'aQ  (52) 


and  its  associated  stretch  (which  we  will  need  for  our  fiber  response  function)  is  still  given  by 

Aa  =  =  A/ao^C'q-o  ■  (53) 

From  equation  52,  we  define  an  intermediate  fiber  a'  in  the  primed  configuration  with  the 
properties  that 

a'  =  F'ao  ,  a  =  Fs  a' .  (54) 

We  can  also  define  stretches  of  a  fiber  in  the  primed  configuration  relative  to  the  other 
configurations.  Let  A'  denote  the  stretch  of  the  primed  fiber  a'  relative  to  the  reference 
configuration  ag: 

A'  =  =  a/ ao^C'ao  ,  (55) 

and  let  As  denote  the  fiber  stretch  of  the  current  configuration  a  relative  to  the  primed 
configuration  a': 


\/ a^a  \/a!'^C sol  Va'^a  ^a^ao  Xa 

\/ a'^ a'  \/ a'^ a'  a/ a'^a'  ^ 


These  are  related  by 


AsA'. 


(56) 


(57) 


Again,  A^  is  the  fiber  stretch  that  results  in  our  final  deformed  configuration  and  that  determines 
the  fiber  response.  A'  is  the  initial  stretch  at  the  start  of  the  simulation,  and  As  is  the  additional 
stretch  in  the  current  configuration  relative  to  the  primed  configuration. 


The  right  Cauchy-Green  deformation  tensor  from  the  simulation,  Cs,  is  calculated  by  SIERRA 
during  the  simulation.  This  means  that  a'  would  need  to  be  determined  at  the  initial  run  time. 
Our  current  method  for  determining  the  fiber  family  orientations  uses  the  primed  coordinate 
system,  so  the  current  algorithm  determines  a  reasonable  a'  direction,  but  with  an  unknown 
magnitude,  since  we  are  requiring  that  |ao|  =  1  and  thus  in  general  |a'|  7^  1.  In  other  words,  if 
we  express  the  Cartesian  coordinates  of  a'  relative  to  a  spherical  coordinate  system. 


a'  =  (A'  cos  9  sin  0,  A'  sin  9  sin  0,  A'  cos  0)  =  X'a  , 


(58) 
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our  algorithm  determines  the  angles  6  and  0  but  A'  is  unknown.  Thus,  if  we  assume  some  F' ,  we 
can  formulate  the  following  constraint  to  determine  A': 


aol  =  \fa^Q  =  ^a'  =  1 . 


(59) 


Then  we  could  use  the  following  expressions  to  calculate  Xa'. 

Xa  =  Vdo^Cao  =  ^Jao^F'^CsF'ao. 


(60) 


A  simpler  alternative  is  to  assume  we  have  an  idea  of  A',  then  we  only  need  to  calculate  A^  to 
determine  A^: 


y/  a^a  \/a''^C  sO,'  a/  X'a^C  sX'a 

\/  a'^a'  \/  a'^  a'  V  X'a^X'a 


(61) 


The  primed  fiber  direction  a  is  already  determined  by  the  pre-existing  algorithm,  and  SIERRA 
provides  Cs,  which  would  play  the  role  of  C  from  section  3.  This  leaves  only  the  user  to  choose 
X'  to  introduce  a  nondimensional  scaling  to  the  fiber  stretch  and  amend  the  implementation  of  the 
fiber  response  function,  equation  30,  F  so  that  A^  is  multiplied  by  the  input  parameter  A',  i.e., 

h  =  {XsX'f. 


6.4  Incorporation  of  Active  Contractile  Fiber  Response 


In  the  previous  section,  we  discussed  a  method  to  introduce  a  prestress  to  the  fiber  families.  In 
this  section,  we  briefly  discuss  how  other  choices  for  the  fiber  response  can  achieve  a  variety  of 
new  behaviors,  such  as  developing  a  simple  constitutive  model  for  skeletal  muscles.  Skeletal 
muscle  has  been  implemented  in  a  variety  of  ways  (29-31).  However,  a  simple  implementation 
that  captures  the  key  features  of  skeletal  muscle  will  be  useful  to  evaluate  whether  such  a  model 
(or  more  complicated  model)  is  needed  in  the  first  place. 

The  fundamental  feature  we  introduce  is  that  in  lieu  of  a  constant  function  -F(Aa),  the  fiber 
response  is  allowed  to  dynamically  evolve.  As  such,  we  present  physical  arguments  that  will 
ultimately  lead  to  an  equation  of  motion  describing  the  dynamic  evolution  of  the  fiber  response 
given  by  F.  One  early  model  for  muscle  was  developed  by  A.  V.  Hill  in  1938  (32).  A 
description  of  Hill’s  model  appears  in  chapter  18  of  the  textbook  by  Keener  and  Sneyd  (33),  the 
key  components  are  only  summarized  here.  The  starting  point  of  this  model  is  the  notion  of  an 
elastic  element  in  series  with  a  contractile  element.  As  such,  the  length  of  a  muscle  fiber  L  is 
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decomposed  into  an  elastic  component  x  and  a  contractile  element  I  so  that  quite  simply, 


L  =  I  +  X  . 


(62) 


By  convention,  the  velocity  of  the  contractile  element  is  given  by  the  negative  time  derivative  of  /, 


V  =  — 


dt 


(63) 


Numerous  experiments  support  a  force- velocity  relationship  between  the  velocity  of  shortening  v 
and  the  load  /: 

if  +  a)  V  =  b{fo-  f)  .  (64) 

Here  a  and  b  are  parameters  to  be  fit  by  the  model  and  /o  is  the  isometric  force,  the  force 
generated  when  the  muscle  length  is  held  fixed.  We  note  that  a  has  dimensions  of  force  and  b  has 
dimensions  of  velocity. 


Assuming  the  force  /  is  a  function  of  the  elastic  length  x,  i.e.,  F{x),  we  can  express  the  time 
derivative  of  /  as 


dF  dx 

dF 

'dL 

dV 

dx  dt 

dx 

dt 

dt 

and  after  inserting  the  force- velocity  relationship,  equation  64,  into  equation  65,  we  obtain 


(65) 


/ 


dF  IdL 

dx  [dt  f  +  a 


(66) 


Following  the  Hill  model,  F  is  taken  to  be  a  linear  function  of  x,  so  that 


F{x)  =  a  {x  —  xq)  , 


(67) 


where  a  is  some  elastic  spring  coefficient, 
muscle  fibers: 


f  =  a 


This  gives  us  the  Hill  model  for  active  contractile 


dL 

dt 


+  b 


fo-f 

f  +  a 


(68) 


At  this  point,  we  introduce  a  possible  extension  of  our  model  to  include  this  behavior.  If,  instead 
of  considering  forces,  we  consider  stresses,  then  we  can  modify  equation  68.  This  modified 
equation  can  then  be  taken  to  directly  affect  the  fiber  response  function  F.  By  letting  f  ^  F, 
a  ^  C,b  ^  /3  and  L  — )■  A^,  we  come  to  the  following  result: 


F 


C 


+  /3 


•Fq-J^ 

F  F  a 


(69) 


30 


C  and  a  have  units  of  stress  and  and  /3  are  strain-rates,  e.g.,  1/s.  This  is  now  an  equation  that 
governs  the  functional  dependence  of  the  fiber  response  that  is  qualitatively  consistent  with  a 
Hill-type  muscle  law.  At  each  time  step,  the  task  is  then  to  calculate  the  muscle  fiber  stretch  \a 
and  its  time  derivative  Aq  so  that  the  fiber  response  at  the  current  time  step  can  be  calculated. 
One  possible  discretization  of  this  is  (first  order,  forward  time): 

^-+1  =J^-  +  C  (Ar^  -  a:)  +  Cp  °  A/ ,  (70) 

y-"-  -f  a 

where  the  superscript  n  denotes  the  time  step  at  which  the  variable  is  evaluated.  While  this  is  not 
a  rigorous  approach  to  incorporating  a  fiber  response  that  mimics  skeletal  muscle,  it  is  an  outline 
of  an  approach  that  could  be  made  rigorous  in  future  studies. 


7.  Concluding  Remarks 


We  have  developed  and  verified  a  transversely  isotropic  hyperelastic  constitutive  model  with  two 
fiber  families  that  can  readily  be  applied  to  soft  biological  tissues  exhibiting  anisotropic  features. 
We  have  outlined  some  current  projects  involving  the  applications  of  this  model,  specifically 
regarding  the  spine.  In  addition,  we  presented  a  clear  pathway  on  how  to  alter  the  preprocessing 
step  so  that  this  model  can  be  used  with  DSI  data  from  brain  tissue  to  capture  a  second  fiber 
direction.  Finally,  we  outlined  a  general  procedure  of  how  the  fiber  response  could  be  easily 
modified  to  represent  prestress  and  active-contractile  muscle  behavior. 
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